MMPPsampler (version 1.0)

Gibbs-sampler for MMPPs: MMPPsampler

Description

Fearnheard and Sherlock (2006) <10.1111/j.1467-9868.2006.00566.x> proposed an exact Gibbs-sampler for performing Bayesian inference on Markov Modulated Poisson processes (MMPP). This package is an efficient implementation of their proposed Gibbs-sampler for binned data generated by an MMPP that uses 'C++' via the 'Rcpp' interface.

Furthermore, the package contains an efficient implementation of the hierarchical MMPP framework, proposed by Clausen (2017) <https://github.com/hc2116/MMPPsampler/blob/master/Master_thesis_Henry.pdf>, that is tailored towards inference on network flow arrival data and extends Fearnheard and Sherlock's Gibbs sampler. This model introduces a latent layer in the arrival process that governs the observed arrival data.

Both sampling frameworks harvests greatly from routines that are optimised for this specific problem in order to remain scalable and efficient for large amounts of input data. These optimised routines include matrix exponentiation and multiplication, and endpoint-conditioned Markov process sampling.

For a detailed description of MMPPS and the exact working of the implemented Gibbs sampler, please see the given references.

Usage

GibbsSampler(y_0T,M,Inter,
             alpha_Gamma_rate, beta_Gamma_rate,
             alpha_Gamma_Q, beta_Gamma_Q,
             alpha_D = NULL, 
             B=20,N=100,
             messages=TRUE)
GibbsSampler_hierarchical(y_0T,M,Inter,
             alpha_Gamma_rate, beta_Gamma_rate,
             alpha_Gamma_Q, beta_Gamma_Q,
             alpha_Gamma_Y, beta_Gamma_Y,
             alpha_D = NULL,
             B=20,N=100,
             messages=TRUE)

Arguments

y_0T

a numeric vector which contains the number of observation in each bin

M

the number of states the unobserved Markov Process \(X_t\) can take

Inter

the temporal length of the binning intervals. All intervals must have the same length. Due to numerical stability, the value for Inter MUST NOT be lower than 1! Please rescale your timescale if values lower are required.

alpha_Gamma_rate

a numeric vector of length \(M\) which specifies the alpha-hyperparameters of the Gamma-prior on the arrival rates \(lambda\).

beta_Gamma_rate

a numeric vector of length \(M\) which specifies the beta-hyperparameters of the Gamma-prior on the arrival rates \(lambda\).

alpha_Gamma_Q

a numeric vector of length \(M\) which specifies the alpha-hyperparameters of the Gamma-prior on the decay rates of \(X_t\) (the diagonal elements of the generator matrix \(Q\).

beta_Gamma_Q

a numeric vector of length \(M\) which specifies the beta-hyperparameters of the Gamma-prior on the decay rates of \(X_t\) (the diagonal elements of the generator matrix \(Q\).

alpha_D

a numeric vector of length \(M-1\) which specifies the alpha-hyperparameters of the Dirichlet-prior on the transition rates of \(X_t\) (the off-diagonal elements of the generator matrix \(Q\). If left blank, will be initialised with 1 (uninformative prior).

alpha_Gamma_Y

a numeric value which specifies the alpha-hyperparameter of the Gamma-prior on the arrival rate \(lambda_Z\) of the latent variable \(Y\).

beta_Gamma_Y

a numeric value of length \(M\) which specifies the beta-hyperparameter of the Gamma-prior on the arrival rate\(lambda_Y\) of the latent variable \(Y\).

B

a numeric value specifying the number of samples discarded as the burn-in of the sampler

N

a numeric value specifying the number of samples to be generated

messages

a boolean value specifying whether to generate text output informing the user about the sampling progress

Value

x

a numeric matrix containing the sampled trajectories of \(X_t\). Each column corresponds to a sample trajectory.

Y

a numeric matrix containing the sampled values of the latent variable \(Y\). Each column corresponds to a sample.

Lambda_S

a numeric matrix containing the sampled values of the arrival rates \(lambda\). Each row corresponds to a sample.

Q_r

a numeric matrix containing the sampled values of the decay rates \(Q_ii\). Each row corresponds to a sample.

Q_d

a numeric matrix containing the sampled values of the transition rates \(Q_i!=j\). Each row corresponds to a sample.

Lambda_Y

a numeric vector containing the sampled values of the arrival rate \(lambda_Y\). Each value corresponds to a sample.

LLH

a numeric vector containing the sampled log-likelihoods of the data with the sampled parameters. Each value corresponds to a sample.

cpu_time

The time needed to create the samples.

Details

Both the Gibbs-sampling functions for the regular MMPP and or the hierarchical extension require an input vector that contains the binned observations, the length of a binning interval, the number of states of the hidden Markov process, and lose prior hyperparameters. As a return, the user receives the desired number of sample trajectories of the hidden Markov process as well as the likelihood of each trajectory.

References

Fearnhead, Paul, and Chris Sherlock. "An exact Gibbs sampler for the Markov-modulated Poisson process." Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68.5 (2006): 767-784 <10.1111/j.1467-9868.2006.00566.x>

Clausen, Henry. "A Bayesian Approach to Human Behaviour Modelling in Computer Networks". Master's thesis, Imperial College London, <https://github.com/hc2116/MMPPsampler/blob/master/Master_thesis_Henry.pdf>

See Also

Source-code and more descriptions available under <https://github.com/hc2116/MMPPsampler>

Examples

# NOT RUN {
#Use the MMPP sample data included in the package to test the regular Gibbs sampler
data("TestdataMMPP")
Test <- TestdataMMPP
test_samples <- GibbsSampler(y_0T=Test$Bins,
                                    M=Test$M,
                                    Inter = Test$Inter,
                                    alpha_Gamma_rate = Test$alpha_Gamma_rate,
                                    alpha_Gamma_Q = Test$alpha_Gamma_Q,
                                    beta_Gamma_Q = Test$beta_Gamma_Q,
                                    beta_Gamma_rate = Test$beta_Gamma_rate,
                                    B=1,N=2,messages=FALSE)

test_plot <- MMPPplot(Sampler_Output = test_samples,
              title = "Example Plot")
plot(test_plot)

#Use the flow sample data included in the package to test the hierarchical model
data("Testdataflows")
Test <- Testdataflows
test_samples <- GibbsSampler_hierarchical(y_0T=Test$Bins,
                                           M=Test$M,
                                           Inter = Test$Inter,
                                           alpha_Gamma_rate = Test$alpha_Gamma_rate,
                                           alpha_Gamma_Q = Test$alpha_Gamma_Q,
                                           beta_Gamma_Q = Test$beta_Gamma_Q,
                                           beta_Gamma_rate = Test$beta_Gamma_rate,
                                           alpha_Gamma_Y=Test$alpha_Gamma_Z,
                                           beta_Gamma_Y=Test$beta_Gamma_Z,
                                           B=1,N=2,messages=FALSE)

# }